#include <cmath>
#include <stdio.h>
#include <iostream>
#include <unistd.h>
#include <time.h>
#include <stdlib.h>
#include "../include/Image.h"
#include "../include/Pixel.h"
#include "../include/vecteur.h"
#include "../include/Masque.h"
#include "../include/Matrice.h"

using namespace std;

int aleat(int min, int max)
{
	return min+(int) (rand() % (max+1-min));
}

Image::Image(int _l,int _h):largeur(_l),hauteur(_h)
{
	data = new Pixel[largeur*hauteur];
	sprite = NULL;
}

Image::~Image(){

	delete []data;
	delete []sprite;
}

Image* Image::CopierImage ( Image * out/*=NULL*/ )
{
	if ( out == NULL )
		out = new Image(largeur,hauteur);
	else if ( out->largeur != largeur || out->hauteur != hauteur )
	{
		std::cerr << "Image::CopierImage () : erreur : images de tailles différentes" << std::endl;
		return this;
	}
	for (int i = 0; i < largeur; ++i){
		for (int j = 0; j < hauteur; ++j){
		   out->data[j*largeur + i] = this->data[j*largeur + i];
		}
	}
	return out;
}

void Image::BlitFrom ( Image * im, int xmin, int xmax, int ymin, int ymax )
{
	if ( im == NULL )
	{
		std::cerr << "Image::BlitFrom () : erreur : pointeur nul donné en paramètre" << std::endl;
		return;
	}
	if ( xmax-xmin < 0 || ymax-ymin < 0 || xmin < 0 || ymin < 0
		|| xmax >= largeur || ymax >= hauteur )
	{
		std::cerr << "Image::BlitFrom () : erreur : rectangle incorrect" << std::endl;
		return;
	}
	for (int y = ymin; y <= ymax; ++y){
		for (int x = xmin; x <= xmax; ++x){
		   this->data[y*largeur + x] = im->data[y*largeur + x];
		}
	}
}

int Image::chargePPM (const char *nom)
{
	FILE *f;
	char s[2048];
	int width, height, maxgrey, n, k;
	unsigned char ligne[2048];
	char * t;
	
	if(data!=NULL)
		delete [] data;
	data = NULL;
	
	/* Ouverture du fichier */
	f = fopen (nom, "rb");
	if (f == NULL)
			{
			printf ("ERREUR fopen(%s)\n", nom);
			return -1;
			}

	/* Lecture MagicValue
	 * On evite fscanf pour eviter buffer overflow */
	t = fgets (s, 2048, f);
	if (s[0] != 'P' || s[1] != '6')
			{
			printf ("ERREUR MagicValue <> P6\n");
			fclose (f);
			return -1;
			}

	/* Lecture width height
	 * On saute les ligne commencant par un '#' */
	do t = fgets (s, 2048, f);
	while (s[0] == '#');
	sscanf (s, "%d %d", &width, &height);
	
	largeur = width;
	hauteur = height;

	/* Lecture maxgrey
	 * On ne tient pas compte de maxgrey, qui est presque toujours 255.
	 * On saute les ligne commencant par un '#' */
	do t = fgets (s, 2048, f);
	while (s[0] == '#');
	sscanf (s, "%d", &maxgrey);
	data = new Pixel[width*height];

	/* Lecture image */
	for (n=0;n<height;n++)
			{
			if(fread(ligne,sizeof(unsigned char)*width*3,1,f)!=1)
				return -1;
			for (k=0;k<width*3;k+=3){
					data[n*width+k/3] = Pixel(ligne[k],ligne[k+1],ligne[k+2]);
				}
			}
	fclose (f);
	return 0;
}

int Image::sauvePPM (const char *nom){
	
	FILE *f;
	int n,k;
	unsigned char ligne[2048];

	f = fopen (nom, "wb");
	if (f == NULL)
	{
		printf ("ERREUR fopen(%s)\n", nom);
		return -1;
	}

	fprintf (f, "P6\n");
	fprintf (f, "# width height\n");
	fprintf (f, "%d %d\n", largeur, hauteur);
	fprintf (f, "255\n");

	for (n=0;n<hauteur;n++)
	{
		for (k=0;k<largeur*3;k+=3){
			ligne[k]=(unsigned char)data[n*largeur+k/3].getR();
			ligne[k+1]=(unsigned char)data[n*largeur+k/3].getV();
			ligne[k+2]=(unsigned char)data[n*largeur+k/3].getB();
		}
		if (fwrite (ligne, sizeof(unsigned char)*largeur*3,1,f)!=1)
		{
			printf ("ERREUR fwrite()\n");
			fclose (f);
			return -1;
		}
	}

	fclose (f);
	return 0;

}

Pixel Image::LitPixel(int x, int y)
{
	if (x < 0)
		x = -x;
	else if (x >= largeur)
		x = largeur*2 - 1 - x;
	if (y < 0)
		y = -y;
	else if (y >= hauteur)
		y = hauteur*2 -1 - y;
	return data[y*largeur + x];
}

void Image::EcritPixel(int x, int y, Pixel couleur)
{
	if (x < 0)
		x = -x;
	else if (x >= largeur)
		x = largeur*2 - 1 - x;
	if (y < 0)
		y = -y;
	else if (y >= hauteur)
		y = hauteur*2 -1 - y;
	data[y*largeur + x] = couleur;
}

Image* Image::FiltrageLineaire(FiltreLineaire *fl){
	Image* out = new Image(largeur,hauteur);
	int bx = (int)fl->taillex/2;
	int by = (int)fl->tailley/2;
	for(int i=0;i<hauteur;i++){
		for(int j=0;j<largeur;j++){
			Pixel s;
			int n=0;
			for(int k=i-by;k<=i+by;k++){
				for(int l=j-bx;l<=j+bx;l++){
					s+=LitPixel(l,k)*fl->coef[n++];
				}
			}
			out->EcritPixel(j,i,s);
		}
	}
	return out;
}

Pixel Image::FiltrageLineaire(FiltreLineaire *fl,int x,int y){
	Pixel res;
	int bx = (int)fl->taillex/2;
	int by = (int)fl->tailley/2;
	int n=0;
	for(int k=y-by;k<=y+by;k++){
		for(int l=x-bx;l<=x+bx;l++){
			res+=LitPixel(l,k)*fl->coef[n++];
		}
	}
	return res;
}


void Image::Afficher(){
	if(sprite!=NULL)
		glDrawPixels(largeur, hauteur, GL_RGB, GL_UNSIGNED_BYTE, sprite) ;
}

void Image::toSprite(){
	
	int x,cpt=0;
	if(sprite != NULL)
		delete []sprite;

	sprite = new GLubyte[hauteur*largeur*3];
	
	for (int y=hauteur-1 ; y > 0; y--)
	{
		for (x=0 ; x < largeur*3 ; x+=3)
		{
			for(int c=0;c<3;c++)
				sprite[cpt++] = (GLubyte)LitPixel(x/3,y)[c] ;

		}
	}
}

Image* Image::Laplacien(){
	Image* laplacien;
	FiltreLineaire fl;
	float tmp[9] = {0,1,0,1,-4,1,0,1,0};
	
	fl.taillex = 3;
	fl.tailley = 3;
	fl.coef = new float[3*3];
	for(int i = 0;i<fl.taillex*fl.tailley;i++){
		fl.coef[i] = tmp[i];
	}
	laplacien = FiltrageLineaire(&fl);
	return laplacien;
}

Pixel Image::Laplacien(int x,int y){

	FiltreLineaire lapl;
	Pixel lapla;
	lapl.taillex = 3;
	lapl.tailley = 3;
	float tmp[9] = {-1,-1,-1,-1,8,-1,-1,-1,-1};

	lapl.coef = new float[9];
	for(int i = 0;i<lapl.taillex*lapl.tailley;i++){
		lapl.coef[i] = tmp[i];
	}

	lapla = FiltrageLineaire(&lapl,x,y);

	delete []lapl.coef;

	return lapla;
}

Image* Image::Gradientx()
{
    Image* gradient;
    FiltreLineaire gradientx;
	float tmp[9] = {1,0,-1,2,0,-2,1,0,-1};

	gradientx.taillex = 3;
	gradientx.tailley = 3;
	gradientx.coef = new float[3*3];
	for(int i = 0;i<gradientx.taillex*gradientx.tailley;i++){
		gradientx.coef[i] = tmp[i];
	}
    gradient = FiltrageLineaire(&gradientx);
    return gradient;
}

Pixel Image::Gradientx(int x,int y){
	int choix = 5;//aleat(0,20);
	Pixel gradient;
	double r,v,b;

	Pixel PC = LitPixel(x-1,y);
	Pixel NC = LitPixel(x+1,y);
	Pixel PP = LitPixel(x-1,y-1);
	Pixel PN = LitPixel(x-1,y+1);
	Pixel NP = LitPixel(x+1,y-1);
	Pixel NN = LitPixel(x+1,y+1);


	if(choix<3){
		/*Sobel etendu*/

		FiltreLineaire gradientx;
		gradientx.taillex = 5;
		gradientx.tailley = 3;
		float tmp[15] = {1,-2,0,2,1,-2,-4,0,4,2,-1,-2,0,2,1};

		gradientx.coef = new float[15];
		for(int i = 0;i<gradientx.taillex*gradientx.tailley;i++){
			gradientx.coef[i] = tmp[i];
		}

		gradient = FiltrageLineaire(&gradientx,x,y);

		delete []gradientx.coef;
	}

	else if(choix >=3 && choix<17){
		/*Central difference*/
		r = (NC[0]-PC[0])/2;
		v = (NC[1]-PC[1])/2;
		b = (NC[2]-PC[2])/2;
		gradient = Pixel(r,v,b);

	}

	else if(choix >=17 && choix<19){
		// Rotation invariant mask.
		const double t = (double)(0.25f*(2-std::sqrt(2.0f))), g = (double)(0.5f*(std::sqrt(2.0f)-1));
		r = -t*PP[0] - g*PC[0] - t*PN[0] + t*NP[0] + g*NC[0] + t*NN[0];
		v = -t*PP[1] - g*PC[1] - t*PN[1] + t*NP[1] + g*NC[1] + t*NN[1];
		b = -t*PP[2] - g*PC[2] - t*PN[2] + t*NP[2] + g*NC[2] + t*NN[2];
		gradient = Pixel(r,v,b);

	}
	else if(choix >=19){
		/*Sobel*/
		r = this->LitPixel(x-1,y-1)[0]-this->LitPixel(x+1,y-1)[0]+2*this->LitPixel(x-1,y)[0]-2*this->LitPixel(x+1,y)[0]+this->LitPixel(x-1,y+1)[0]-this->LitPixel(x+1,y+1)[0];
		v = this->LitPixel(x-1,y-1)[1]-this->LitPixel(x+1,y-1)[1]+2*this->LitPixel(x-1,y)[1]-2*this->LitPixel(x+1,y)[1]+this->LitPixel(x-1,y+1)[1]-this->LitPixel(x+1,y+1)[1];
		b = this->LitPixel(x-1,y-1)[2]-this->LitPixel(x+1,y-1)[2]+2*this->LitPixel(x-1,y)[2]-2*this->LitPixel(x+1,y)[2]+this->LitPixel(x-1,y+1)[2]-this->LitPixel(x+1,y+1)[2];
		gradient = Pixel(r,v,b);
	}
	return gradient;
}

Image* Image::Gradienty()
{
    Image* gradient;
    FiltreLineaire gradienty;
	float tmp[9] = {-1,-2,-1,0,0,0,1,2,1};

	gradienty.taillex = 3;
	gradienty.tailley = 3;
	gradienty.coef = new float[3*3];
	for(int i = 0;i<gradienty.taillex*gradienty.tailley;i++){
		gradienty.coef[i] = tmp[i];
	}

    gradient = FiltrageLineaire(&gradienty);
    return gradient;
}

Pixel Image::Gradienty(int x,int y){

	int choix = 5;//aleat(0,20);
	Pixel gradient;

	double r,v,b;

	Pixel PP = LitPixel(x-1,y-1);
	Pixel PN = LitPixel(x-1,y+1);
	Pixel NP = LitPixel(x+1,y-1);
	Pixel NN = LitPixel(x+1,y+1);
	Pixel CN = LitPixel(x,y+1);
	Pixel CP = LitPixel(x,y-1);

	if(choix<3){

		FiltreLineaire gradienty;
		gradienty.taillex = 3;
		gradienty.tailley = 5;
		float tmp[15] = {-1,-2,-1,-2,-4,-2,0,0,0,2,4,2,1,2,1};

		gradienty.coef = new float[15];
		for(int i = 0;i<gradienty.taillex*gradienty.tailley;i++){
			gradienty.coef[i] = tmp[i];
		}

		gradient = FiltrageLineaire(&gradienty,x,y);

		delete []gradienty.coef;
	}

	else if(choix >=3 && choix<17){
		/*Central difference*/
		r  = (CN[0]-CP[0])/2;
		v= (CN[1]-CP[1])/2;
		b= (CN[2]-CP[2])/2;
		gradient = Pixel(r,v,b);
	}

	else if(choix >=17 && choix<19){
		// Rotation invariant mask.
		const double t = (double)(0.25f*(2-std::sqrt(2.0f))), g = (double)(0.5f*(std::sqrt(2.0f)-1));
		r = -t*PP[0] - g*CP[0] - t*NP[0] + t*PN[0] + g*CN[0] + t*NN[0];
		v = -t*PP[1] - g*CP[1] - t*NP[1] + t*PN[1] + g*CN[1] + t*NN[1];
		b = -t*PP[2] - g*CP[2] - t*NP[2] + t*PN[2] + g*CN[2] + t*NN[2];
		gradient = Pixel(r,v,b);
	}
	else if(choix >=19){
		/*Sobel*/
		r = -this->LitPixel(x-1,y-1)[0]-2*this->LitPixel(x,y-1)[0]-this->LitPixel(x+1,y-1)[0]+this->LitPixel(x-1,y+1)[0]+2*this->LitPixel(x,y+1)[0]+this->LitPixel(x+1,y+1)[0];
		v = -this->LitPixel(x-1,y-1)[1]-2*this->LitPixel(x,y-1)[1]-this->LitPixel(x+1,y-1)[1]+this->LitPixel(x-1,y+1)[1]+2*this->LitPixel(x,y+1)[1]+this->LitPixel(x+1,y+1)[1];
		b = -this->LitPixel(x-1,y-1)[2]-2*this->LitPixel(x,y-1)[2]-this->LitPixel(x+1,y-1)[2]+this->LitPixel(x-1,y+1)[2]+2*this->LitPixel(x,y+1)[2]+this->LitPixel(x+1,y+1)[2];
		gradient = Pixel(r,v,b);
	}
	return gradient;
}

vecteur* Image::ComputeSmoothness( int x, int y )
{
	vecteur* smoothness;
	smoothness = new vecteur[3];

	Pixel xm1 = this->Laplacien(x+1,y);
	Pixel xp1 = this->Laplacien(x-1,y);
	Pixel ym1 = this->Laplacien(x,y+1);
	Pixel yp1 = this->Laplacien(x,y-1);

	smoothness[0].x = (xp1[0] - xm1[0]);
	smoothness[0].y = (yp1[0] - ym1[0]);

	smoothness[1].x = (xp1[1] - xm1[1]);
	smoothness[1].y = (yp1[1] - ym1[1]);

	smoothness[2].x = (xp1[2] - xm1[2]);
	smoothness[2].y = (yp1[2] - ym1[2]);

	smoothness[0].normalise();
	smoothness[1].normalise();
	smoothness[2].normalise();
	return smoothness;
}

vecteur* Image::computeIsophoteDirection(int _x,int _y,int inv){
	vecteur* N = new vecteur[3];


	N[0] = vecteur(-Gradienty(_x,_y)[0],Gradientx(_x,_y)[0]);
	N[1] = vecteur(-Gradienty(_x,_y)[1],Gradientx(_x,_y)[1]);
	N[2] = vecteur(-Gradienty(_x,_y)[2],Gradientx(_x,_y)[2]);

	N[0].normalise();
	N[1].normalise();
	N[2].normalise();


	return N;
}

Image* Image::DiffusionAnisotrope (Masque* m)
{
	Image* I=NULL;
	Image* I1=NULL;
	Image* swap;
	I = new Image(this->largeur,this->hauteur);
	int lambda = 30;
	double dt = 0.25;
	double gradNR, gradSR, gradER, gradWR;
	double gradNV, gradSV, gradEV, gradWV;
	double gradNB, gradSB, gradEB, gradWB;
	double cNR, cSR, cER, cWR;
	double cNV, cSV, cEV, cWV;
	double cNB, cSB, cEB, cWB;
	double r,v,b;

	I1 = this->CopierImage();
	I = this->CopierImage();

	list<Point>::iterator it;

	for ( int nb = 0 ; nb < 5 ; nb++ )
	{
		swap = I1;
		I1 = I;
		I=swap;
		for(it = m->listePos.begin();it!=m->listePos.end();it++){

			int i = (*it).x;
			int j = hauteur -(*it).y;
			/* Calcul des 4 gradients */
			gradNR = I->LitPixel(i,j-1).getR() - I->LitPixel(i,j).getR();
			gradSR = I->LitPixel(i,j+1).getR() - I->LitPixel(i,j).getR();
			gradER = I->LitPixel(i+1,j).getR() - I->LitPixel(i,j).getR();
			gradWR = I->LitPixel(i-1,j).getR() - I->LitPixel(i,j).getR();

			gradNB = I->LitPixel(i,j-1).getB() - I->LitPixel(i,j).getB();
			gradSB = I->LitPixel(i,j+1).getB() - I->LitPixel(i,j).getB();
			gradEB = I->LitPixel(i+1,j).getB() - I->LitPixel(i,j).getB();
			gradWB = I->LitPixel(i-1,j).getB() - I->LitPixel(i,j).getB();

			gradNV = I->LitPixel(i,j-1).getV() - I->LitPixel(i,j).getV();
			gradSV = I->LitPixel(i,j+1).getV() - I->LitPixel(i,j).getV();
			gradEV = I->LitPixel(i+1,j).getV() - I->LitPixel(i,j).getV();
			gradWV = I->LitPixel(i-1,j).getV() - I->LitPixel(i,j).getV();

			/* calcul des 4 coefficients de diffusivite */
			cNR = exp ( - ( (gradNR*gradNR) / (lambda*lambda) ) );
			cSR = exp ( - ( (gradSR*gradSR) / (lambda*lambda) ) );
			cER = exp ( - ( (gradER*gradER) / (lambda*lambda) ) );
			cWR = exp ( - ( (gradWR*gradWR) / (lambda*lambda) ) );

			cNV = exp ( - ( (gradNV*gradNV) / (lambda*lambda) ) );
			cSV = exp ( - ( (gradSV*gradSV) / (lambda*lambda) ) );
			cEV = exp ( - ( (gradEV*gradEV) / (lambda*lambda) ) );
			cWV = exp ( - ( (gradWV*gradWV) / (lambda*lambda) ) );

			cNB = exp ( - ( (gradNB*gradNB) / (lambda*lambda) ) );
			cSB = exp ( - ( (gradSB*gradSB) / (lambda*lambda) ) );
			cEB = exp ( - ( (gradEB*gradEB) / (lambda*lambda) ) );
			cWB = exp ( - ( (gradWB*gradWB) / (lambda*lambda) ) );

			r = I->LitPixel(i,j).getR() + dt * ( cNR * gradNR + cSR * gradSR + cER * gradER + cWR * gradWR );
			v = I->LitPixel(i,j).getV() + dt * ( cNV * gradNV + cSV * gradSV + cEV * gradEV + cWV * gradWV );
			b = I->LitPixel(i,j).getB() + dt * ( cNB * gradNB + cSB * gradSB + cEB * gradEB + cWB * gradWB );

			I1->EcritPixel(i,j, Pixel (r,v,b) );
		}
	}

	delete I;

	return I1;
}

Matrice* Image::computeTenseurStructure ( int _x, int _y, int _canal )
{
	Matrice* resultat = new Matrice(2,2);
	Matrice* gradient = new Matrice(1,2);
	Matrice* transpose ;
	Pixel tmp;
	
	if ( _canal == 0 )
	{
		/* on travaille sur le rouge */
		
		/* on remplit le gradient */
		tmp = this->Gradientx(_x,_y);
		gradient->matrice[0][0] = tmp.getR();
		tmp = this->Gradienty(_x,_y);
		gradient->matrice[0][1] = tmp.getR();
	}
	
	if ( _canal == 1 )
	{
		/* on travaille sur le rouge */
		
		/* on remplit le gradient */
		tmp = this->Gradientx(_x,_y);
		gradient->matrice[0][0] = tmp.getV();
		tmp = this->Gradienty(_x,_y);
		gradient->matrice[0][1] = tmp.getV();
	}
	
	if ( _canal == 2 )
	{
		/* on travaille sur le rouge */
		
		/* on remplit le gradient */
		tmp = this->Gradientx(_x,_y);
		gradient->matrice[0][0] = tmp.getB();
		tmp = this->Gradienty(_x,_y);
		gradient->matrice[0][1] = tmp.getB();
	}
	
	
	/* on calcule sa transposée */
	transpose = gradient->transposee(gradient);
	
	resultat = gradient->multiplication(transpose);
	return resultat;
}

Matrice* Image::computeTenseurDiffusion ( int _x, int _y, int _canal )
{
	Matrice* tenseurStructure;
	Matrice* vecteurPropre;
	Matrice* matricePropre; // matrice obtenue par multiplication du vecteur propre et de sa transposée
	Matrice* resultat;
	float coeff;
	
	tenseurStructure = this->computeTenseurStructure(_x,_y,_canal);
	vecteurPropre = tenseurStructure->vecteurPropreMin();
	matricePropre = vecteurPropre->multiplication(vecteurPropre->transposee());
	
	coeff = 1.0/( sqrt(1+tenseurStructure->valeurPropreMax()+tenseurStructure->valeurPropreMin() ) );
	
	resultat = matricePropre->multiplication( coeff );
	
	return resultat;
}
/*
Matrice * Image::computeHessienne ( int _x, int _y, int _canal/ *= 0* /,
	Matrice * _mat/ *= NULL* / )
{
	if ( _canal < 0 || _canal >= 3 )
		_canal = 0;
	if ( _mat == NULL )
		_mat = new Matrice ( 2, 2 );
	Image * _img = this;
	
	_mat->setValeur ( 0, 0,
		_img->LitPixel ( _x+1, _y )[_canal]
		+ _img->LitPixel ( _x-1, _y )[_canal]
		- 2 * _img->LitPixel ( _x, _y )[_canal] );
	_mat->setValeur ( 1, 0,
		0.25 * ( _img->LitPixel ( _x+1, _y+1 )[_canal]
		+ _img->LitPixel ( _x-1, _y-1 )[_canal]
		- _img->LitPixel ( _x-1, _y+1 )[_canal]
		- _img->LitPixel ( _x+1, _y-1 )[_canal] ) );
	_mat->setValeur ( 0, 1, _mat->getValeur ( 1, 0 ) );
	_mat->setValeur ( 1, 1,
		_img->LitPixel ( _x, _y+1 )[_canal]
		+ _img->LitPixel ( _x, _y-1 )[_canal]
		- 2 * _img->LitPixel ( _x, _y )[_canal] );
	return _mat;
}
*/

Matrice * Image::computeHessienne ( int _x, int _y, int _canal/*= 0*/,
	Matrice * _mat/*= NULL*/ )
{ // fonction très sollicitée, doit être la plus rapide possible
	if ( _canal < 0 || _canal >= 3 )
		_canal = 0;
	if ( _mat == NULL )
		_mat = new Matrice ( 2, 2 );
	
	// pour éviter de passer par LitPixel () qui fait perdre du temps
	if ( _x < 1 )
		_x = 1;
	else if ( _x > largeur - 2 )
		_x = largeur - 2;
	if ( _y < 1 )
		_y = 1;
	else if ( _y > hauteur - 2 )
		_y = hauteur - 2;
	
	double v0 = (data[(_y-1)*largeur + (_x-1)])[_canal];
	double v1 = (data[(_y-1)*largeur + (_x)])[_canal];
	double v2 = (data[(_y-1)*largeur + (_x+1)])[_canal];
	double v3 = (data[(_y)*largeur + (_x-1)])[_canal];
	double v4 = (data[(_y)*largeur + (_x)])[_canal];
	double v5 = (data[(_y)*largeur + (_x+1)])[_canal];
	double v6 = (data[(_y+1)*largeur + (_x-1)])[_canal];
	double v7 = (data[(_y+1)*largeur + (_x)])[_canal];
	double v8 = (data[(_y+1)*largeur + (_x+1)])[_canal];
	
	_mat->setValeur ( 0, 0,
		v5
		+ v3
		- 2.0 * v4 );
	_mat->setValeur ( 1, 0,
		0.25 * ( v8
		+ v0
		- v6
		- v2 ) );
	_mat->setValeur ( 0, 1, _mat->getValeur ( 1, 0 ) );
	_mat->setValeur ( 1, 1,
		v7
		+ v1
		- 2.0 * v4 );
	return _mat;
}

Pixel Image::operator()(const int x, const int y){
	return LitPixel(x,y);
}

Matrice * Image::computeGradient ( int _x, int _y, int _canal/*= 0*/,
	Matrice * _mat/*= NULL*/ )
{ // fonction très sollicitée, doit être la plus rapide possible
	if ( _canal < 0 || _canal >= 3 )
		_canal = 0;
	if ( _mat == NULL )
		_mat = new Matrice ( 1, 2 );
	
	// pour éviter de passer par LitPixel () qui fait perdre du temps
	if ( _x < 1 )
		_x = 1;
	else if ( _x > largeur - 2 )
		_x = largeur - 2;
	if ( _y < 1 )
		_y = 1;
	else if ( _y > hauteur - 2 )
		_y = hauteur - 2;
	
	double v0 = (data[(_y-1)*largeur + (_x-1)])[_canal];
	double v1 = (data[(_y-1)*largeur + (_x)])[_canal];
	double v2 = (data[(_y-1)*largeur + (_x+1)])[_canal];
	double v3 = (data[(_y)*largeur + (_x-1)])[_canal];
	//double v4 = (data[(_y)*largeur + (_x)])[_canal];
	double v5 = (data[(_y)*largeur + (_x+1)])[_canal];
	double v6 = (data[(_y+1)*largeur + (_x-1)])[_canal];
	double v7 = (data[(_y+1)*largeur + (_x)])[_canal];
	double v8 = (data[(_y+1)*largeur + (_x+1)])[_canal];
	
	_mat->setValeur ( 0, 0,
		- v0
		- 2 * v3
		- v6
		+ v2
		+ 2 * v5
		+ v8
		 );
	_mat->setValeur ( 0, 1,
		- v0
		- 2 * v1
		- v2
		+ v6
		+ 2 * v7
		+ v8
		 );
	return _mat;
}

Matrice * Image::nablaInablaIT ( int _x, int _y, int _canal/*= 0*/, Matrice * _res/*= NULL*/ )
{
	if ( _canal < 0 || _canal >= 3 )
		_canal = 0;
	Image * _img = this;
	Matrice * gradient = _img->computeGradient ( _x, _y, _canal, NULL );
	_res = gradient->mmT2 ( _res ); // _res sera créé s'il était NULL
	delete gradient;
	return _res;
}

Matrice * Image::sommeNablaInablaIT ( int _x, int _y, Matrice * _m22temp, Matrice * _res/*= NULL*/,
	int _canaux/*= 3*/ )
{
	/* on supprime les tests de validité pour gagner du temps
	if ( _canaux <= 0 || _canaux > 3 )
		_canaux = 1;
	if ( _m22temp == NULL )
	{
		std::cerr << "Image::sommeNablaInablaIT () : erreur : pointeur nul donné en paramètre" << std::endl;
		return NULL;
	}*/
	//Matrice * ggTCur = new Matrice ( 2, 2 );
	if ( _res == NULL )
		_res = new Matrice ( 2, 2 );
	_res->setValeur ( 0, 0, 0.0 );
	_res->setValeur ( 0, 1, 0.0 );
	_res->setValeur ( 1, 0, 0.0 );
	_res->setValeur ( 1, 1, 0.0 );
	for ( int i = 0; i < _canaux; ++i )
	{
		//nablaInablaIT ( _x, _y, i, ggTCur );
		nablaInablaIT ( _x, _y, i, _m22temp );
		//_res->addition ( ggTCur, _res );
		_res->addition ( _m22temp, _res );
	}
	//delete ggTCur;
	return _res;
}
